%% Plotting setting
%clear;
set(0, 'DefaultAxesFontSize',30);
set(0, 'DefaultLineLineWidth',6);
%% Psycometric and chronometric function
%load('monk_data.mat');

%Experimental Monkey data
c0 = [0.032 0.064 0.128  0.256 0.512]';
RT1_EXP = [795.0 819.2; 758.5 837.7; 671.8 776.7; 532.4 539.0; 409.3 NaN];
RT1_STD = [ 12.8 15.8; 10.4 19.0; 8.3 43.4; 6.1 0; 4.9 NaN];
Prop1 = [ 337 193; 409 125; 486 32; 511 1;512 0];
RT0_EXP = [ 818.3 872.5; 758.3 823.7; 678.1 890.7; 551.1 785.3; 436.9 NaN];
RT0_STD = [ 12.5 18.4; 10.7 23.4; 8.8 42.4; 6.1 21.1; 4.6 NaN];
Prop0 = [323 175; 387 104; 477 28; 510 4; 516 0];




% c0 = task(1).coh_set;
% Prop0 = task(1).pT1;
% RT0 = [task(1).rtc; task(1).rtw]';
% RT0_STD = [task(1).rtc_se; task(1).rtw_se]';
%load decision Boundary
d0 = load('policy_-0.1_50.0_0.0.txt');

nT = round(size(d0, 1)* 0.8);

   
[~, RT0] = modelRT(d0(1:nT,:),c0);
%Find out scale for reaction time and non-decision time
linSqrFit =  [RT0(:,1), ones(length(c0),1)] \ RT0_EXP(:,1);

c = [exp(log(0.01):0.2:log(0.512)) 0.512]';
[PR, RT0] = modelRT(d0(1:nT,:),c);

 
RT0 = RT0 * linSqrFit(1) + linSqrFit(2);

figure;
subplot(1,2,1);
plot(c, PR);
hold on;
plot(c0, Prop0(:,1)./ (Prop0(:,1)+Prop0(:,2)),'ok',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',10);
xlabel('Coherence','FontWeight', 'bold');
ylabel('Proportion Correct','FontWeight', 'bold');
xlim([0.01, 0.512]);
ylim([0.5 1]);
set(gca, 'XScale', 'log', 'XTick',[0.01 0.1 0.5], 'LineWidth',2,...
    'FontWeight','bold');
subplot(1,2,2);
plot(c, RT0(:,1), 'b');
hold on;
errorbar(c0, RT0_EXP(:,1), RT0_STD(:,1),'ok',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',10);
%errorbar(c0, RT0(:,2), RT0_STD(:,2),'or',...
%     'MarkerFaceColor',[0 0 0],'MarkerSize',10);
%legend('correct choices \pm SEM', 'Incorrect choices \pm SEM',...
%   'Location','SouthWest');
xlabel('Coherence', 'FontWeight', 'bold','FontSize',30);
ylabel('Reaction Time (ms)', 'FontWeight', 'bold','FontSize',30);
xlim([0.01, 0.512]);
set(gca, 'XScale', 'log', 'XTick',[0.01 0.1 0.5], 'LineWidth',2,...
    'FontWeight','bold');
set(gcf,'paperunits','inches');
set(gcf,'papersize',[18 12]);
set(gcf,'paperposition',[0,0,18,12]);
% saveas(gcf,'PCRT.fig','fig');
% saveas(gcf,'PCRT.eps','epsc2');

%% PCRTs
d3 = load('policy_-0.1_1000.0_0.0.txt');
d2 = load('policy_-0.1_100.0_0.0.txt');  
d1 = load('policy_-0.1_50.0_0.0.txt');
d0 = load('policy_-0.1_25.0_0.0.txt');

nT = round(size(d0, 1)* 0.8);

c = -0.5:0.01:0.5;
[PR0, RT0] = modelRT(d0(1:nT,:),c);

[PR1, RT1] = modelRT(d1(1:nT,:),c);
[PR2, RT2] = modelRT(d2(1:nT,:),c);
[PR3, RT3] = modelRT(d3(1:nT,:),c);

RT0(c < 0,:) = fliplr(RT0(c < 0,:));
RT1(c < 0,:) = fliplr(RT1(c < 0,:));
RT2(c < 0,:) = fliplr(RT2(c < 0,:));
RT3(c < 0,:) = fliplr(RT3(c < 0,:));
 
figure;
subplot(1,2,1);
plot(c, PR0, '-','DisplayName', 'Ratio = 250');
hold on;
plot(c, PR1, '-g','DisplayName', 'Ratio = 500');
plot(c, PR2, '-r','DisplayName', 'Ratio = 1000');
plot(c, PR3, '-c','DisplayName', 'Ratio = 10000');
hold off;
xlabel('Coherence','FontWeight', 'bold');
ylabel('Proportion Correct','FontWeight', 'bold');
xlim([0 0.512]); 
ylim([0.5 1]);
legend_handle = legend('show','Location', 'Best');
set(legend_handle, 'Box', 'off')
set(gca, 'LineWidth',2, 'FontWeight','bold');

subplot(1,2,2);
plot(c, RT0(:,1),'-','DisplayName', 'Ratio = 250');
hold on;
plot(c, RT1(:,1), '-g','DisplayName', 'Ratio = 500');
plot(c, RT2(:,1), '-r','DisplayName', 'Ratio = 1000');
plot(c, RT3(:,1), '-c','DisplayName', 'Ratio = 10000');
xlabel('Coherence', 'FontWeight', 'bold','FontSize',30);
ylabel('POMDP Steps', 'FontWeight', 'bold','FontSize',30);
xlim([-0.5, 0.5]);
ylim([0,320]);
legend_handle = legend('show', 'Location', 'NorthWest');
set(legend_handle, 'Box', 'off');
set(gca,  'LineWidth',2,...
    'FontWeight','bold');
% 
% 
set(gcf,'paperunits','inches');
set(gcf,'papersize',[18 12]);
set(gcf,'paperposition',[0,0,18,12]);
% saveas(gcf,'PCRTs.fig','fig');
% saveas(gcf,'PCRTs.eps','epsc2');



%%  Speed and Accuracy Trade-off
 
ssp = load('Lsp_data.mat');
sacc = load('Lacc_data.mat');

d1 = load('policy_-0.1_125.0_0.0.txt');
d0 = load('policy_-0.1_45.0_0.0.txt');

nT = round(size(d0, 1)* 0.8);

c = ssp.task.coh_set;

[PR0, RT0] = modelRT(d0(1:nT,:),c);

[PR1, RT1] = modelRT(d1(1:nT,:),c);

RT0(c < 0,:) = fliplr(RT0(c < 0,:));
RT1(c < 0,:) = fliplr(RT1(c < 0,:));


linSqrFit =  [[RT0(:,1);RT1(:,1)], ones(length(c) * 2,1)] ...
    \ [ssp.task(1).rtc'; sacc.task(1).rtc'];

linSqrFit0 =  [RT0(:,1), ones(length(c),1)] \ ssp.task(1).rtc';
linSqrFit1 =  [RT1(:,1), ones(length(c),1)] \ sacc.task(1).rtc';


RT0 = RT0 * linSqrFit(1) + linSqrFit(2);

RT1 = RT1 * linSqrFit(1) + linSqrFit(2);

 
figure;
subplot(1,2,1);
plot(c, PR0);
hold on;
plot(c, ssp.task(1).pT1,'ok',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',25);
plot(c, sacc.task(1).pT1,'xr',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',25);
plot(c, PR1,'--r');
xlabel('Coherence','FontWeight', 'bold');
ylabel('Proportion Correct','FontWeight', 'bold');
xlim([0 0.512]); 
ylim([0.5 1]);
set(gca, 'LineWidth',2,'FontWeight','bold');

subplot(1,2,2);
plot(c, RT0(:,1));
hold on;
plot(c, RT1(:,1),'--r');
errorbar(c, ssp.task(1).rtc, ssp.task(1).rtc_se,'ok',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',25);
errorbar(c, sacc.task(1).rtc, sacc.task(1).rtc_se, 'xr',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',25);
xlabel('Coherence', 'FontWeight', 'bold','FontSize',30);
ylabel('Reaction Time (ms)', 'FontWeight', 'bold','FontSize',30);
xlim([-0.5, 0.5]); 
ylim([300, 2000]);
set(gca, 'YScale', 'log', 'YTick',[400 600 800 1000 1500], 'LineWidth',2,...
    'FontWeight','bold');
% 
% 
% set(gcf,'paperunits','inches');
% set(gcf,'papersize',[18 12]);
% set(gcf,'paperposition',[0,0,18,12]);
% saveas(gcf,'SpeedAccuracy.fig','fig');
% saveas(gcf,'SpeedAccuracy.eps','epsc2');

%% Effects of Prior Probability on reaction times

load('monk_data.mat');
c = task(1).coh_set;

d0_B = load('policy_biased.txt');
%plotPolicy(d0_B);
% saveas(gcf,'policy_biased.fig','fig');
% saveas(gcf,'policy_biased.jpg','jpg');

d0_U = load('policy_-0.1_50.0_0.0.txt');


[PR_U, RT_U] = modelRT(d0_U,c);
[PR_B, RT_B] = modelRT(d0_B,c);

RT_U(c < 0, :) = fliplr(RT_U(c<0,:));
RT_B(c < 0, :) = fliplr(RT_B(c<0,:));

linSqrFit_U =  [RT_U(:,1), ones(length(c),1)] \ task(1).rtc';

linSqrFit_B =  [RT_B(:,1), ones(length(c),1)] \ task(2).rtc';

RT_U = RT_U * linSqrFit_U(1) + linSqrFit_U(2);
RT_B = RT_B * linSqrFit_B(1) + linSqrFit_B(2);
%RT_B = RT_B * linSqrFit_U(1) + linSqrFit_U(2);

figure;
subplot(2,1,1);
plot(c, PR_U, '-');
hold on;
plot(c, PR_B, '--g');
plot(c, task(1).pT1,'o',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',10);
plot(c, task(2).pT1,'og',...
    'MarkerFaceColor',[0 0 0],'MarkerSize',10);
hold off;
xlim([min(c) max(c)]);
ylabel('Proportional Correct','FontWeight', 'bold','FontSize',30);
set(gca, 'XTick', [-.5 -.25 0 .25 .5],'XMinorTick', 'on',  'LineWidth',2,...
    'FontWeight','bold')
subplot(2,1,2);
%Show reaction times of correct trials only
plot(c, RT_U(:,1),'-');
hold on;
plot(c, RT_B(:,1),'--g');
errorbar(c, task(1).rtc, task(2).rtc_se, 'o', ...
    'MarkerFaceColor',[0 0 0],'MarkerSize',10);
errorbar(c, task(2).rtc, task(2).rtc_se, 'go', ...
    'MarkerFaceColor',[0 0 0],'MarkerSize',10);
hold off;
xlim([min(c) max(c)]);
xlabel('Motion Strength','FontWeight', 'bold','FontSize',30);
ylabel('Reaction Time','FontWeight', 'bold','FontSize',30);
set(gca, 'XTick', [-.5 -.25 0 .25 .5],'XMinorTick', 'on', 'LineWidth',2,...
    'FontWeight','bold');
set(gcf,'paperunits','inches');
set(gcf,'papersize',[12 15]);
set(gcf,'paperposition',[0,0,12,15]);
% saveas(gcf,'prior_monk.fig','fig');
% saveas(gcf,'prior_monk.jpg','jpg');